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Abstract 

Presented here is an algorithm for a type-II quantum computer which simulates the 
Ising model in one and two dimensions. It is equivalent to the Metropolis Monte- 
Carlo method and takes advantage of quantum superposition for random number 
generation. This algorithm does not require the ensemble of states to be measured 
at the end of each iteration, as is required for other type-II algorithms. Only the 
binary result is measured at each node which means this algorithm could be im- 
plemented using a range of different quantum computing architectures. The Ising 
model provides an example of how cellular automata rules can be formulated to be 
run on a type-II quantum computer. 
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1 Introduction 



There has been much effort recently in the development of quantum com- 
puting schemes for a range of applications. One such scheme consists of an 
array of small quantum information processors connected by classical com- 
munication channels and is termed a type-II quantum computer (T2QC) [1]. 
Figure 1 shows conceptually the design of a T2QC where each node contains a 
small number of qubits (10's or less) and the nodes are connected via classical 
communication channels. Each node is initialised into some quantum state 
and then the entire computer undergoes some global unitary operation. The 
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state of each node is then measured and the results used to re-initialise the 
neighbouring nodes for the next computational step. These are known as the 
collision and streaming steps respectively, using the lattice-gas terminology. 
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Fig. 1. Design of a Type-II Quantum Computer, composed of many small conven- 
tional quantum computers 

This form of quantum computer is able to solve certain fluid dynamics prob- 
lems which are difficult to solve on a classical computer [2] [3] [4]. While this 
is a useful application in its own right, it also provides an important 'step- 
ping stone' towards globally phase coherent quantum computers of the type 
required to run Shor's factoring algorithm [5] and other important quantum 
computing algorithms. 

The type-II architecture is also well suited to a range of other physical prob- 
lems. One such system is the Metropolis Monte-Carlo (MC) method for sim- 
ulating the Ising model. It has been previously suggested that this problem 
can be simulated using a set of cellular automaton rules on a lattice [6] . Using 
this approach, the problem is well suited to a type-II quantum computer as 
it requires many lattice points but only nearest neighbour interactions. While 
there are several other deterministic or semi-deterministic schemes available to 
simulate the Ising model with simple cellular automata [7] [8] [9], they are only 
approximations to the metropolis solution. This is because the random num- 
ber generation required by the metropolis algorithm is difficult to implement 
without float-point numbers or at least very large integers. A cellular automa- 
ton rule typically requires only small numbers of bits or integers. The major 
advantage of using a T2QC to implement the metropolis algorithm is that the 
random number generation can be included using quantum superposition. 

In this paper, the required evolution rules are developed for the Metropolis 
algorithm and then the necessary unitary evolution is determined using the 
quantum circuit formalism. The probabilistic aspect of the this algorithm is 
incorporated in the form of the weighted coin-toss application of quantum su- 
perposition. The necessary quantum evolution circuits are derived for both the 
one-dimensional (ID) and two-dimensional (2D) Ising models. The implemen- 
tation of these circuits is discussed for technologies where only the resulting 
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binary state can be measured and also where the ensemble of states can be 
measured. 



2 Metropolis simulation of the Ising model 



The Ising model consists of a regular array of spins Sj which take one of two 
values {—1,1}. The energy at each site % is given by 

Ei = -J^SiSj (1) 



where the sum is over nearest neighbours (j — i±l) and the total energy of the 
lattice is given by the sum of the on-site energies. The use of periodic boundary 
conditions is assumed, though fixed boundaries can be easily implemented. As 
the spins can take only two values (assuming there is no external field) the 
energy at each site can only take a finite number of values, the number of which 
depend on the dimension of the lattice. For the case of a ID chain of Ising 
spins, the change in energy can take three distinct values, AE = OJ, ±4J. 
The value of J is assumed to be constant across the lattice and for most of the 
discussion, it is assumed to be positive (J > 0). This situation corresponds 
to the ferromagnetic case, though the generalisation to the anti-ferromagnetic 
case (J < 0) is relatively simple. 

Metropolis Monte-Carlo [10] is a popular method of simulating the Ising model 
for a given temperature. The metropolis algorithm involves flipping a spin 
at random and calculating the corresponding change in energy {AE). This 
spin flip is then accepted with probability given by min{l,e~~kr}, where T is 
the temperature of the lattice and k is the Boltzmann constant. Throughout 
this discussion it is assumed that k = 1 which results in temperature being 
expressed in units of J. 

As the energy levels given by equation 1 are discrete, the metropolis algo- 
rithm is easily converted to a set of cellular automaton rules. In order to 
implement this, the same evolution rule must be applied to the entire lat- 
tice simultaneously. However, if the metropolis algorithm is applied to every 
site simultaneously the 'feedback catastrophe' results, as pointed out by Vich- 
niac [6]. Instead, the rule must be applied to (at most) every second site in 
a checkerboard configuration [11]. This can be achieved by either storing two 
sublattices or using a parity bit which controls whether the evolution rule is 
applied or not [7]. 
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3 Quantum circuit formalism 



As the quantum circuit formalism is universal and encompasses classical boolean 
logic [12], type-II quantum computers can be designed to simulate most de- 
terministic and probabilistic cellular automata. The probabilistic aspect can 
be incorporated through the use of a qubit in a superposition state to provide 
a weighted coin-toss probability. The Ising model algorithm is presented as an 
example of this idea, though the process should generalise for most cellular 
automata rules. 

The appropriate evolution rules required to implement the Metropolis algo- 
rithm are given here using the conventional quantum circuit notation [12] [13]. 
This formalism provides a simple way to describe the evolution while offering a 
straight forward method of checking the operation via a truth-table approach. 
The notation used to represent the reversible gates used in these circuits is 
briefly described here. For a more complete discussion, see chapters 1 and 4 
of Nielsen and Chuang[13]. 

The CNOT or controlled-NOT gate is illustrated in figure 2 and is accom- 
panied by the corresponding truth-table. The basic operation of the CNOT 
gate is to invert the target bit (\b)) if the control bit (\a)) is equal to one. The 
convention is to use a filled circle to represent the control bit while inversion 
of the target bit is indicated by the 'crossed-circle'. The CNOT gate can be 
extended to include two or more input states, an example of which is the CC- 
NOT or Toffoli gate, shown in figure 3. Similarly the Toffoli gate inverts the 
target bit if both control bits are equal to 1. Other gates of interest are the 
inversion or NOT gate (figure 4(a)) and the swap gate (figure 4(b)) which just 
swaps the input states \a) and \b). If the dependance on the control bit of any 
of these gates is reversed, this is indicated by a open circle as shown in figure 
5. The fact that the control bits for both the CNOT and Toffoli gates have 
not changed after the operation is a result of the reversibility of these gates. 

While these reversible gates are essentially a classical concept, the power of 
quantum computing arises when they are applied to quantum superposition 
states where a qubit can be in some superposition of and 1. In order to 
prepare these superposition states, the quantum computer must be able to 
perform so called 'rotations' to set up an arbitrary superposition. The result 
of this is that the control and target qubits for the quantum version of these 
gates can be in some unknown state which is only determined by measurement 
of the qubit. The measurement process returns one of two possibilities where 
the probability of returning a given state is the square of the amplitude of 
that state in the superposition. The application of this kind of measurement is 
indicated by the symbol in figure 4(c) and is termed a destructive measurement 
as all information about the quantum superposition is lost. 
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Fig. 2. Notation for a CNOT gate and the truth-table showing classical operation 
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Fig. 3. Notation for a Toffoli gate and the truth-table showing classical operation 



\a) — — \a) 



(a) 



|a>— ¥ 



|6> 
|a> 



(b) 



(c) 



Fig. 4. Notation for a NOT gate (a), a swap gate (b) and a destructive measurement 
(c) 
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Fig. 5. Notation for a CNOT gate where the target bit is inverted if the control bit 
is equal to 



For each circuit in this paper, \S)i is designated the on-site Ising spin (sj) where 
Si = and \S')i is the new on-site spin after application of the quantum 
circuit. This allows a simple correspondence between the Ising states and the 
qubit encoding, giving 'spin-up' (+1) as binary state 1 and 'spin-down' (— 1) 
as state 0. 
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4 ID Ising model 



For a particular on-site spins Si of the ID Ising model, the neighbouring spins 
are designated A and B. As the spins only take discrete values the change in 
energy due to a single spin flip can only take a finite set of values (AE = OJ, 
±4 J), as long as there is no global field. This means the Metropolis method 
can be implemented using the circuit given in figure 6, where \P) is a qubit 
initialised in the state \P) = y/P\l) + y/l - P|0) with P = e" giving the 
appropriate probability for the given temperature and coupling strength. The 
circuit also has an 'ancilla' qubit which is initialised with the state |0) and is 
used for temporary storage during the computation. 

The general algorithm steps are as follows 

(1) Initialise on-site spin for each node (\S)i) based on the initial state of the 
Ising lattice 

(2) Initialise input states A and B based on the value of the neighbouring 
spins using the streaming rule (see section 6) 

(3) Apply quantum circuit gate sequence to the entire computer simultane- 
ously 

(4) Measure resultant spin for each node 

(5) Use resultant spin to re-initialise neighbouring nodes 

(6) Repeat until equilibrium is reached 
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Fig. 6. Quantum circuit for the evolution of the on-site spin \S) for the ID Ising 
model 

The parameters of interest (mean magnetisation for instance) can be calcu- 
lated once the output of each node is measured. In figure 6 the circuit requires 
5 qubits and 4 multi-qubit gates though it is possible to reduce the num- 
ber of qubits by rewriting the circuit in terms of two-qubit gates and single 
qubit rot at ions [12]. The choice of which two-qubit gates to use for this de- 
composition is largely a matter of which gates are most suitable for a given 
implementation. In this form the circuit will require at least 2d+ 1 qubits for a 
ci-dimensional lattice to represent the on-site spin and its nearest neighbours. 

The truth-table given in figure 7 shows that this circuit is equivalent, once 
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measured, to the classical Metropolis method. This table gives the output of 
the quantum circuit \S') for each input state ASB as well as the classical result 
S' cl and its probability of acceptance p c i. The result of the circuit given in 6, 
once the state \S') is measured, is identical to the classical result. All but two 
of the entries in this table are exactly equivalent to the classical metropolis 
result, where a spin flip is accepted as long as the change in energy is zero or 
negative. The first (last) entry differs in that the result is a superposition state 
with a probability P of obtaining a spin-up (spin-down) result and a corre- 
sponding probability 1 — P of obtaining the opposite result. This is equivalent 
to the random number generation step of the metropolis algorithm where a 
configuration is rejected or accepted based a weighted probability. 
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Fig. 7. Truth-table for the quantum circuit in 6, showing the equivalence to the 
classical Metropolis algorithm once the superposition state |S") is measured 

5 2D Ising model 

Using the same approach, the Metropolis method for the 2D Ising model can 
be implemented using the circuit given in figure 8. 

Inputs A-D are the neighbouring spins and Pi and P 2 are the probabilities 
P 1 = e~^r~ and P 2 = e~~ respectively. These correspond to qubits initialised 
in the states |Pi) = y/l\\l) + VI - A|0) and \P 2 ) = ^/P~ 2 \l) + VI - P 2 |0). 

The circuit given for the 2D Ising model is not optimal and can in fact be 
realised with considerably fewer and simpler gates using a 'full-adder' type 
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Fig. 8. Quantum circuit for 2D Ising model 

configuration. It is presented in this form as it is more apparent which sections 
of the circuit operate for AE = J, ±4 J and ±8 J, making the operation of the 
circuit more obvious. If this circuit were to be realised experimentally, it could 
be broken down into the appropriate two qubit operations which suited the 
particular implementation [12]. This can result in considerable simplification, 
depending on the details of the experimental system. 



6 Streaming and parallelisation 



The streaming process for the ID Ising chain is illustrated in figure 9. Each 
node is initialised using the value for the on-site spin Si and copies of its 
nearest neighbours Ai = Si-i and Bi = Si+i. The superposition qubit (P) 
and the ancilla qubit (An) for each node are re-initialised every iteration with 
the same value, as discussed in section 4 and 5. 




Node i-1 Node i Node i + 1 



Fig. 9. Each node is initialised with the on-site spin Si and copies of the neighbouring 
on-site spins S*j_i and Si + i 

The streaming rule for the 2D Ising lattice is similar to the ID case. As dis- 
cussed earlier, all nodes cannot be updated simultaneously otherwise the lat- 
tice will oscillate and not reach a stable state. A checkerboard update scheme 
can be used instead, as illustrated in figure 10. Each 'black' node is initialised 
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with its spin value Sij and its four nearest neighbours A = Sij+i, B = Si-ij, 
C = Sij-i and D = S^+ij. The new value of S is calculated and then the 
same process is followed for the 'white' nodes. Once this process is completed, 
every node in the lattice has been updated. 
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Si+lj 
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Fig. 10. The checkerboard update scheme with the spin value Sij and its four nearest 
neighbours A = <Sij+i, B = -S^-ij, C = -S'jj-i and D = Si+ij illustrated 

As a result of the reversibility of quantum circuits and the fact that the tar- 
get qubits are left unchanged, the values of the neighbouring spins A-D are 
preserved after the computation. This means that no additional storage space 
is required to store the spins that are not being modified at each evolution 
(assuming error free operation). It also means that the value of an arbitrary 
neighbour (ie: A) can be used to initialise the spin S at the next time step. 
This results in only half as many nodes being required to simulate the lattice. 
This algorithm therefore requires y computational nodes for N Ising spins. 
In this situation two evolutions of the quantum circuit are required to update 
every node in the lattice once. An alternative solution is to include a parity 
bit which controls whether the rest of the circuit functions or not. The lattice 
is then initialised so that the parity bits alternate across the lattice and they 
are inverted at each evolution. 



7 Possible implementations 

The algorithm, as discussed, is independent of implementation and would 
therefore be suited to a range of different quantum computing architectures. 
Of particular interest is that the final state of the spin after each evolution 
collapses to either a 1 or a when measured. Existing diffusion based type-II 
algorithms require the probability amplitudes to be measured at each step, 
which requires either time or ensemble averaging. This cellular automaton 
based algorithm, on the other hand, does not require the measurement of 
the full ensemble of states and is therefore suited to solid state or ion-trap 
architectures. 

One important consideration is that the accuracy with which the temperature 
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of the system can be controlled is directly related to the accuracy with which 
the superposition state \P) can be created. The amplitude (p = \fP) of this 
superposition is given by p = e~ 2 ^ T for the 2D case and p = e~~ x l T for ID. The 
accuracy with which this superposition can be created and controlled puts a 
lower bound on the temperature resolution the system can simulate. This is 
especially important for architectures where only a binary output is measured, 
as the superposition must be created using single qubit rotations. The accuracy 
in the temperature is therefore related to the accuracy of the rotation gate. By 
using standard uncertainty analysis, the maximum allowable error rate (dp) 
can be calculated for a given temperature resolution using equation 2. 



Sp 



dp 
dT 



ST 



(2) 



Figure 11 shows the maximum allowed gate error rate for a single qubit rota- 
tion (Sp) needed to achieve a temperature resolution (ST) of 0.1, as a function 
of system temperature (T). The curves p\ and p 2 correspond to the maxi- 
mum allowed gate error rate (Sp) for the ID and 2D systems respectively. For 
T < 0.5 the required gate accuracy is effectively infinite but in this region 
the simple Ising system displays trivial behaviour. For the region of interest 
(1.5 < T < 4), a maximum error rate of Sp = 0.01 is sufficient for ID and 2D 
calculations with a temperature resolution of 0.1. The two-qubit gates must 
also be controlled with a comparable resolution, though the net effect of these 
gates should be at most linear in the number of gates. 
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Fig. 11. Maximum allowable error rate required for a temperature resolution of 
5T = 0.1, as a function of Ising lattice temperature 
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8 Ensemble streaming 



The algorithm detailed so far is formally equivalent to the classical metropolis 
algorithm for the Ising model. An interesting side note involves determining 
what would happen if the measurement step were to be replaced with an en- 
semble measurement. An ensemble measurement consists of measuring many 
copies of the system and determining the magnitude of the probability distri- 
bution for each qubit. This distribution can then be used to re-initialise the 
other spins, ready for the next iteration of the algorithm. This situation cor- 
responds to either time averaging over many runs or ensemble averaging, as 
is the case with NMR based implementations, and will be referred to as 'en- 
semble streaming'. This is especially interesting as the only type-II quantum 
computer in existence at the time of writing is based on the NMR architec- 
ture and therefore uses ensemble streaming[14] [15]. The result of this process 
is that the spins no longer take on definite integer values but can have con- 
tinuous values depending on their neighbours and the temperature dependant 
function. The magnetisation as a function of temperature is shown for this 
ensemble streaming (Ensemble T2QC) for the 2D Ising model (the dot-dash 
line in figure 12). 

By allowing the full ensemble of states to propagate through the lattice, the 
model now bears some similarity to a mean-field style approximation. As each 
node can vary continuously from to 1, the deviation from the mean spin 
value quickly decays to zero for any one spin site. After a few iterations, 
each node on the lattice takes on the same expectation value, resulting in 
a smooth curve even for a very small lattice of four spins. As a result, the 
curve for any number of spins is found to be the same as for the four spin 
case. In addition the ensemble version converges much faster than the classical 
Metropolis algorithm (less than 5 iterations per temperature point). The curve 
(figure 12) corresponds to the system being heated from the ground state up 
to a state where every configuration is equally likely, which is equivalent to 
the randomised state. If the system was then cooled again, it tends to stay 
completely randomised and therefore the zero magnetisation state is the stable 
solution all the way back to zero temperature. This is a consequence of allowing 
the spins to be continuous (an ensemble of states) rather than taking discrete 
values and is a non-physical solution. 

For comparison, the solution for the ensemble T2QC is plotted in figure 12 
with the analytic solution due to Onsager [16] for an infinite lattice, a clas- 
sical Metropolis MC run using a 1000 x 1000 point lattice and the classical 
mean field solution for an infinite lattice [17]. The classical Metropolis MC 
run was performed using a temperature step size of 0.01 with a minimum of 
20 and a maximum of 10000 iterations per temperature point. This includes 
1000 2 spin flips per iteration with each spin flip requiring the generation of a 
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random number, the calculation of the change in energy and the evaluation 
of an exponential weighting function. By comparison, for the standard T2QC 
algorithm with one-shot measurement all these steps are performed through 
the evolution of the circuit in figure 8. This occurs for each spin in the lattice 
simultaneously which means that the circuit need only be applied once per it- 
eration. The number of iterations required for convergence using the one-shot 
measurement algorithm is similar to the classical algorithm. 
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Fig. 12. Magnetisation as a function of temperature for the 2D Ising model. The 'Ex- 
act' solution is the analytic solution for any infinite lattice, 'Mean field' is the mean 
field approximation for an infinite lattice, 'MC is a classical Metropolis monte-carlo 
run using a 1000 x 1000 lattice and 'Ensemble T2QC is a simulated run of the Ising 
algorithm for a T2QC using ensemble streaming on 4 nodes 



The ensemble measurement algorithm underestimates the critical tempera- 
ture, as can be seen in figure 12, giving T c m 2.109, compared to the analytic 
solution of T c = 2.269 [16]. The form of the ensemble solution is a polynomial 
approximation rather than the accepted power law behaviour. This is due to 
the fact that as the true Ising lattice approaches the critical temperature, 
large clusters of spins form in random locations in the lattice, resulting in 
non-zero magnetisation. The ensemble measurement algorithm, on the other 
hand, effectively averages over all possible configurations. As the clusters oc- 
cur in random locations their effects cancel out, resulting in a zero average 
magnetisation. This discrepancy accounts for the difference in critical temper- 
ature and gives some insight into the regime where clustering is an important 
effect. It is expected that a large cluster generalisation of this algorithm would 
improve the critical behaviour near T c as the forming of clusters would be more 
accurately modelled. 
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9 Conclusion 

A novel algorithm for a type-II quantum computer has been presented for sim- 
ulating the Metropolis Monte-Carlo method for the Ising model in one and 
two dimensions. The algorithm requires only N/2 nodes to simulate N spins 
and is formally equivalent to a probabilistic cellular automaton formulation of 
the Metropolis method for the Ising model. It requires a maximum of 5 qubits 
per node for the ID case and 12 qubits per node for the 2D case, though this 
is an upper bound for arbitrary interactions. This compares well with other 
classical implementations which require many more classical bits for random 
number generation. It also demonstrates that other cellular automata models 
may be simply reformulated to run on a type-II quantum computer. This ex- 
tends the range of applications for a type-II quantum computer beyond those 
based on the lattice Boltzmann hydrodynamic equations, such as diffusion. 
Additionally, the effect of running this algorithm on a type-II quantum com- 
puter of the type demonstrated experimentally, using ensemble streaming, was 
studied. While this ensemble streaming variation of the algorithm was found 
to replicate the phase change of a two dimensional Ising lattice, it is unable 
to accurately describe the effects of clustering near this transition. 
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